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Abstract 

We study a fractional reaction-diffusion system with two types of variables: acti- 
vator and inhibitor. The interactions between components are modeled by cubical 
nonlinearity. Linearization of the system around the homogeneous state provides 
information about the stability of the solutions which is quite different from linear 
stability analysis of the regular system with integer derivatives. It is shown that by 
combining the fractional derivatives index with the ratio of characteristic times, it 
is possible to find the marginal value of the index where the oscillatory instability 
arises. The increase of the value of fractional derivative index leads to the time 
periodic solutions. The domains of existing of periodic solutions for different pa- 
rameters of the problem are obtained. A computer simulation of the corresponding 
nonlinear fractional ordinary differential equations is presented. It is established by 
computer simulation that there exists a set of stable spatio-temporal structures of 
the one-dimensional system under the Neumann and periodic boundary condition. 
The features of these solutions consist in the transformation of the steady state 
dissipative structures to homogeneous oscillations or space temporary structures at 
a certain value of fractional index and the ratio of characteristic times of system. 

Key words: reaction-diffusion system, fractional differential equations, 
oscillations, dissipative structures, pattern formation, spatio-temporal structures 



1 Introduction 

Reaction-diffusion systems (RDS) have been extensively used in the study of 
self-organization phenomena in physics, biology, chemistry, ecology etc. (See, 
for example [1,2,3,4, 5,6,7,8,9, 10J). The main result obtained from these sys- 
tems is that nonlinear phenomena include diversity of stationary and spatio- 
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temporary dissipative patterns, oscillations, different types of waves, excitabil- 
ity, bistability etc. The mechanism of the formation of such type of nonlinear 
phenomena and the conditions of their emergence have been extensively stud- 
ied during the last couple decades. 

In the recent years, there has been a great deal of interest in fractional reaction- 
diffusion (FRD) systems [TT1I12|I13II14II15|I16II17II18|I19| which from one side ex- 
hibit self-organization phenomena and from the other side introduce a new 
parameter to these systems, which is a fractional derivative index, and it gives 
a great degree of freedom for diversity of self-organization phenomena. 

At the same time, analyzing structures in FRD systems evolves, both from the 
point of view of the qualitative analysis and from the computer simulation. 
Namely these two problems are the goals of our investigation here. Our partic- 
ular interest is the analysis of the specific non-linear system of FRD equations. 
We consider a very well-known example of the RDS with cubical nonlinearity 
[314TT] which probably is the simplest one used in RD systems modeling 

Ti ^M = pvVM) +m - nl/3 - n 2 , (1) 



subject to 
(i) Neumann: 



T J^n 2 &t) = L 2 V 2 M ^ t) _ n2 + /3ni + A (2) 



drii/dx\ x=0 = drii/ dx\ x= i x = 0,i = 1,2. (3) 



or 



(ii) Periodic: 

rii(t, 0) = rii(t,l x ), drii/dx\ x= Q = dni/dx\ x= i x ,i = 1,2. (4) 

boundary conditions and with the certain initial condition ni\ t= o = n^(x). Here 
x : < x < l x ; (x, t) G R x R + ; V 2 = J^;ni(x,t), n 2 (x,t) G R - activator 
and inhibitor variables correspondingly, n,T2,l,L, G R - characteristic times 
and lengths of the system, A G R - is an external parameter. 

Fractional derivatives 9 q^'^ on the left hand side of equations ([I]),© instead 
of standard time derivatives are the Caputo fractional derivatives in time of 
the order < a < 2 and are represented as [2~Tf2"2] 



T 



8 a 1 r n< m \ 

dt^ 71 ^ : ~ r(m-a) J (t - r)"+ 1 - m 



dr, m — 1 < a < m,m E N. 
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It should be noted that equations ([I]),© at a = 1 correspond to standard 
RDS. At a < 1, they describe anomalous sub-diffusion and at a > 1 
anomalous super-diffusion. 



2 Linear stability analysis 



Stability of the steady-state constant solutions of the system ([I]),© corre- 
spond to homogeneous equilibrium state 

W = m-nl/3-n 2 = 0, Q = -n 2 + fim +A = (5) 

can be analyzed by linearization of the system nearby this solution. In this 
case the system ([I]),© can be transformed to linear system at equilibrium 
point §5§ and linearized around this equilibrium state. As a result we have: 

; =F(u)u(x,t), (6) 



dt a 

/ 



(Z 2 V 2 + aii)/n a 12 /n 

a 2 i/r 2 (L 2 V 2 + a 2 2)/r 2 



Ani(x,t) 

where u(x, t) = \ | , F{u) 

An 2 (x,t) 

is a Frechet derivative with respect to u(x,t), an = W' nx , a\ 2 = W^ 2 , a 2 \ = 
Q' ni i a 22 = Qn 2 ( au derivatives are taken at homogeneous equilibrium states 

( Ani(t)\ 

([5])). By substituting the solution in the form u(x, t) = \ COS rCX, rv — 

\An 2 (t)) 

f-j,j = 1, 2, ... into FRD system (J6j) we can get the system of linear ordinary 
differential equations ([6]) with the matrix F determined by the operator F. 

For analyzing stability conditions of the equations ([I]),© let us use simple 
linear transformation which can convert this linear system ([6]) to a diagonal 
form 



( Ai \ 

where C is a diagonal matrix for F: C = P FP = , eigenvalues Ai 2 

\ A 2 ) 

are determined by the characteristic equation of the matrix F, Ai )2 = |(trF± 

I Ani(t) 

ytr 2 F — 4detF), r/(t) = P \ , P is the matrix of eigenvectors of 



An 2 (t) 



matrix F. 



3 



\ 


a ' 


detF 

/ 






} 


■ 




\ 


1 






\ \ 


/ 






\ \ 


/ 






\ 1 


i / 






\ 


\ / 






\ 








\ 
















Im X 



jarg(X) 



trF 



ReX 



(a) (b) 

Fig. 1. Schematic view of the marginal curve of a (solid line), describing fixed points 
for two-dimensional vector field - (a), the position of eigenvalue A corresponding to 
marginal value of a in the coordinate system (ReX,ImX) - (b). Shaded domains 
correspond to instability region 



In this case, the solution of the vector equation 
functions [20ll2T1l22ll23ll24l 



is given by Mittag-Leffler 



oo /\ £Ot\k 

An *(*) = £ tvV I ^ A ^(°) = E a (X t t a )A ni {0),t = 1,2. (8) 
^T(ka + 1) 

Using the result obtained in the papers [TrJinZE] . we can conclude that if for 
any of the roots 

\Arg(Xi)\ <mr/2 (9) 

the solution has an increasing function component then the system is asymp- 
totically unstable. 

Analyzing the roots of the characteristic equations, we can see that at 4 det F— 
tr 2 F > eigenvalues are complex and can be represented as 



Ai, 2 = -(trF ± iV4 detF-ir 2 F) = \ Re ± i\ Im . 

The roots Ai^ are complex inside the parabola (Fig. Ufa)) and the fixed points 
are the spiral source (trF > 0) or spiral sinks (trF < 0). The plot of the 
marginal value a : a = = -\Arg(Xi)\ which follows from the conditions (Q 
is given by the formula 

I ^ arctan A /4det F/tr 2 F - 1, trF > 0, 

<*>={* V , " (10) 

I 2 - I arctan ^4 det F/tr 2 F - 1, trF < 

and is presented in the Fig. [D 

Let us analyze the system solution with the help of the Fig.[T](a). Consider the 
parameters which keep the system inside the parabola. It is a well-known fact, 
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that at a = 1 the domain on the righthand side of the parabola (trF > 0) 
is unstable with the existing limit circle, while the domain on the left hand 
side (trF < 0) is stable. By crossing the axis trF = the Hopf bifurcation 
conditions become true. 

In the general case of a : < a < 2 for every point inside the parabola 
there exists a marginal value of ao where the system changes its stability. The 
value of a is a certain bifurcation parameter which switches the stable and 
unstable state of the system. At lower a : a < ao = -\Arg(Xi)\ the system 
has oscillatory modes but they are stable. Increasing the value of a > ao = 
-\Arg(\i) \ leads to instability. As a result, the domain below the curve a , as 
a function of trF is stable and the domain above the curve is unstable. 

The plot of the roots, describing the mechanism of the system instability, can 
be understood from the Fig. QJb) where the case «o > 1 is described. In fact, 
having complex number A^ with ReXi < at a — > 2 it is always possible 
to satisfy the condition \Arg(Xi)\ < «7r/2, and the system becomes unstable 
according to homogeneous oscillations (Figdfb)). The smaller is the value of 
trF, the easier it is to fulfill the instability conditions. 

In contrast to this case, a complex values of A,, with ReXi > lead to the sys- 
tem instability for regular system with a — 1. However fractional derivatives 
with a < 1 can stabilize the system if a < cto = -\Arg(\i)\. This makes it 
possible to conclude that fractional derivative equations with a < 1 are more 
stable that their integer twinges. 



3 Solutions of the coupled fractional ordinary differential equations 
(FODEs) 

Let us first consider the coupled fractional ordinary differential equations 
(FODEs) which can be obtained from (JT|),(j2DatZ = L = and analyze the 
stability conditions for such systems. The plot of isoclines for this model is 
represented on Fig. [2](a). In this case homogeneous solution can be determined 
from the system of equations W = Q = and is given by the solution of cubic 
algebraic equation 

(P - l)nx + n?/3 + A = 0. (11) 
Simple calculation makes it possible to write useful expressions required for 
our analysis 



-/3/r 2 1/t 2 J r i r 2 nr 2 
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Fig. 2. Null isoclines - (a), 3-d instability domains in coordinates (ao, n\,T\lr-i) 
- (b), dependence of t\/t2 on n\ at «o changing from 0.1 (bottom curve) to 1.0 
(upper curve) with step 0.1 - (c), dependence of T\/t2 on n\ at ao changing from 
1.1 (bottom curve) to 1.9 (upper curve) with step 0.1 - (d). The shaded domains 
correspond to those one obtained by slicing 3-d surface represented on figure (b). 

It is easy to see that if the value of t\/t^ in certain cases, is smaller than 1, the 
instability conditions (trF > 0) lead to Hopf bifurcation for regular system 
(a = 1) [T|2|3|lf5] . In this case, the plot of the domain, where instability exists, 
is shown on the Fig. [2](c) (unstable domain is below the upper curve). The 
linear analysis of the system for a = 1 shows that, if T\/t2 > 1, the solution 
corresponding to the intersections of two isoclines is stable. The smaller is 
the ratio of t\/t2, the wider is the instability region. Formally, at t\/t2 — > 0, 
the instability region in fix coincides with the interval (—1, 1) where the null 
isocline W(n\,n2) = has its increasing part. These results are very widely 
known in the theory of nonlinear dynamical systems [Tp|3Pf5] . 

In the FODEs the conditions of the instability change (Q, and we have to 
analyze the real and the imaginary part of the existing complex eigenvalues, 
especially the equation: 

4 det F - tr 2 F = 4((/3 - 1) + nj)/^ - ({1 -n\) /r x - l/r 2 f > 0. (12) 

In fact, with the complex eigenvalues, it is possible to find out the correspond- 
ing value of a where the condition ([9]) is true. We will show that this interval 
is not correlated with the increasing part of the null isocline of the system. 
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Indeed, omitting simple calculation, we can write an equation for marginal 
values of H\ 



n\ - 2(1 + -)n\ + -\ - 2— (20 - 1) + 1 = 0, 




(13) 



and solution of this quadratic equation gives the domain where the oscillatory 
instability arises: 



This expression estimates the maximum and minimum values of H\ where the 
system can be unstable at certain value of a = a® as a function of ti/t%. On 
the Fig. 2 the dependance Txjri is given as a function of n\ for different values 
of a changing with the step 0.1. 

In fact, examine the domain of the FODEs where the eigenvalues are complex. 
The condition flUD determine dependence among values rTi,Ti/V 2 and ao (this 
dependance is represented on Fig. 2(b). On the Fig. 2(c) and 2(d) cross sec- 
tions of this figure is presented for fixed value of olq. Inside the curve system of 
FODEs is unstable and outside it is stable. We can see that at a < instabil- 
ity domain in coordinates ( n\, t\/t<i) is smaller then for the case olq = 1 (Fig. 
2(c)). In the same time for the case a > instability domain increases and 
the greater is the value of t\/ti then the wider is the region in n\ where the 
instability holds true (Fig. 2(d)). In this case, from (fill) we immediately ob- 



tain approximate dependance (ni ~ ±v/ti/t 2 ). Of course, at r 1 /r 2 ,n 1 ^> 1 the 
interval in n x where the instability arises increases but the domain decreases 
and gets the form of a narrow concave stripe following parabola ti/t 2 — nf 
for 2 - a < 1. 

Finally, we can conclude that for t\/t2 > 1 in contrast to the regular system 
with integer index a fractional differential equations can be both stable and 
unstable. The greater is the value of ti/t 2 then the wider is interval in n\ where 
instability conditions are true. Opposite situation is for a < 1 while regular 
system with integer index is always unstable fractional dynamical system can 
be stable. It is a statement, that FODEs are at least as stable as their integer 
order counterparts |22j. It is obvious that this statement is true for a < 1 
only. 

It should be noted that, even if the eigenvalues are not complex (Aj m = 0), the 
systems with fractional derivatives can poses oscillatory damping oscillations. 
Such situation takes place when 4detF — tr 2 F < 0, trF < 0, det F > 
and two eigenvalues are real and less than zero. In this case, at 1 < a < 2 
steady state solutions of the system are stable and any perturbations are 
damping. Such system was considered, for example, in the article [25], where 
an analytical solution for fractional oscillator is obtained. 




(14) 
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Our task here is to confirm our linear analysis by finding out not only the 
conditions of the bifurcation but also the real time dynamics of FODEs by 
corresponding combination of the parameters T\/r2,n\ for different values of 
a : < a < 2. We have established here that the dynamics of the FODEs can 
be much more complicated than that one of the equations with integer order 
[2611271128] . 

Fig. 3 gives the results of computer simulation of fractional ODEs for different 
cases considered above. We single out two different cases a < 1 and a > 1. 
The characteristic instability domain (shaded region) for these two cases lo- 
cated on the first row in two columns. Outside the shaded region system is 
stable and solution is given by the equations ([SD- For any point inside shaded 
region system is unstable and we investigate nonlinear dynamics by computer 
simulation of the FODEs. We present nonlinear dynamics for several points 
denoted by capital letters A,B,C,D,E corresponding to homogeneous distribu- 
tions in these points (null isoclines intersect in these points). For the points A 
and B taken from the certain domain obtained at a = 0.8 and for points C,D,E 
taken from the instability domain obtained at a=1.8 nonlinear dynamics is 
presented on Fig. 3(b)-(f). 

From Fig. 3, we can tell that at the point A increment of the oscillations is 
small enough and formation of the limit cycle is during long period (t ~ 100). 
At the point B oscillations develop more rapidly and have a bigger amplitude. 
The result of computer simulation of the system with T\jri > 1 is represented 
in points C, D, E. In this case in point C oscillations have a small amplitude 
and sufficiently large transient time. In the points D, E oscillations get fast 
enough and have a great amplitude. 

We presented simulation of the evolutionary dynamics for only two values of 
indexes a. Nevertheless, similar dynamics is inherent to any value of a (we 
investigated the alpha from a=0.1 till a— 1.9 with step 0.1). 



4 Computer simulation of pattern formation 



This section contains a discussion of the results of the numerical study of the 
initial value problem of the system ([I]),©. The system with corresponding 
initial and boundary conditions was integrated numerically using the explicit 
and implicit schemes with respect to time and centered difference approxi- 
mation for spatial derivatives. The fractional derivatives were approximated 
using the scheme on the basis of Grunwald-Letnikov definition for < a < 2 
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Fig. 3. Instability domains for a = 0.8 and a = 1.8 - (a). The time domain oscil- 
lations (left) and corresponding two dimensional phase portrait (right) for point A 
- (b), B - (c), C - (d), D - (e) E - (f) (r 2 = 1,/? = 2,1 = L = 0). For each point 
A, B, C, D and E parameter A corresponds to steady state solution taken in these 
points. 
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[22] . In other words, for the system of n fractional RD equations 



T J dta] - = d 3 Q X 2 + fj{u 1 ,...,U n ), j = l,Tl, (15) 

where Tj, dj, fj - certain parameters and nonlinearities of the RD system cor- 
respondingly. In this case Grunwald-Letnikov scheme can be represented as 

k _ djjAtp , k k k x {At) a i f ( k k 

U j,i Tj(Ax) 2 W'*- 1 3,i + U hi+y T . J3\ U l,ii ■■■> U n,i) — 



k 



c^ } = l, ^ = ^(1-^), 1 = 1,2, 



where = Uj(xi, t k ) = Uj(iAx, kAt), m 



a. 



The applied numerical schemes are implicit, and for each time layer they 
are presented as the system of algebraic equations solved by Newton- Raphson 
technique. Such approach makes it possible to get the system of equations with 
band Jacobian for each node and to use the sweep method for the solution 
of linear algebraic equations. Calculating the values of the spatial derivatives 
and corresponding nonlinear terms on the previous layer, we obtained explicit 
schemes for integration. Despite the fact that these algorithms are quite sim- 
ple, they are very sensitive and require small steps of integration, and they 
often do not allow to find numerical results. In contrast, the implicit schemes, 
in certain sense, are similar to the implicit Euler's method, and they have 
shown very good behavior at the modeling of fractional reaction-diffusion sys- 
tems for different step size of integration, as well as for nonlinear function and 
the power function of fractional index. 

We have considered here the kinetics of formation of dissipative structures for 
different values of a. These results are presented on Fig. H]and[5j 

The simulations were carried out for a one-dimensional system on an equidis- 
tant grid with spatial step h changing from 0.01 to 0.1 and time step At 
changing from 0.001 to 0.1. We used imposed Neuman (J3]) or periodic bound- 
ary conditions (jl]). As the initial condition, we used the uniform state which 
was superposed with a small spatially inhomogeneous perturbation. 

It is very well known that in the case of a = 1 at l/L 1, trF < 0, det F < 
the homogeneous distribution is unstable according to wave number 



k = k = (aua 2 2 - ai 2 a 2 i) ' (LI) 



l/4/r n -l/2 
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Fig. 4. Numerical solution of the fractional reaction-diffusion equations (|T|) , (J2J) . 
Dynamics of variable m (left column) and n-2 (right column) on the time inter- 
val (0,20) for a = 0.8, k = 8, A = -0.1, fi = 2,r 2 = 1, I 2 = 0.05, L 2 = 1; 
t\/t2 = 0.75 - (a), t\/t2 = 0.37 - (b), t\/t2 = 0.3 - (c). Initial conditions are: 

Til = Hi — 0-05cOs(kQx), rig = U2- 

As a result at certain parameter A such that Hi G (—1,1) the system can 
be unstable and the steady state solutions in the form of nonhomogeneous 
dissipative structures arise. Such type systems have rich dynamics, including 
steady state dissipative structures, homogeneous and nonhomogeneous oscil- 
lations, and spatiotemporal patterns. In this paper, we focus mainly on the 
study of general properties of the solutions depending on the value of a and 
the ratio of the characteristic times t\/t2- 

Fig. H] (a)-(c) show evolution of the dissipative structure formation for a < 1. 
At certain values of T\/ti we have steady state solution - (a), then with de- 
creasing ratio of Ti/t 2 nonhomogeneous pulsating structures - (b) and trans- 
formation of this oscillating regime to homogeneous oscillations - (c). 

The evolutionary dynamics for a > 1 as the value of T\jri decreases is shown 
on the Fig. [3(a)- (c). First plot corresponds to the steady state structures - (a) 
which start to oscillate at smaller ratio of characteristic times and eventually 
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Fig. 5. Numerical solution of the fractional reaction-diffusion equations (|T|) , (J2J) . 
Dynamics of variable ri\ (left column) and n 2 (right column) on the time inter- 
val (0,80) for a = 1.8, k = 4,A = -0.1, f3 = 2,r 2 = l,l 2 = 0.05, L 2 = 1; 
t\/t2 = 10 - (a), t\/t2 = 7 - (b), T1/T2 = 6 - (c). Initial conditions are: 
n\ = ni — 0.05cos(kox),ri2 = n 2 . 

homogeneous oscillation are taking place in FRD system - (c). Such behav- 
ior is due to the case that the oscillatory perturbations are damping in the 
first situation, and then small oscillations are steady state. With the further 
decreasing of the value ti/t 2 the steady state oscillations are unstable and 
dynamics changes to the homogenous temporary behavior (Fig. [3(c)). 

The emergence of homogeneous oscillations, which destroy pattern formation 
(Fig. E^a)-(c)) has deep physical meaning. The matter is that the station- 
ary dissipative structures consist of smooth and sharp regions of variable ni, 
and the smooth shape of n 2 . The linear system analysis shows that the ho- 
mogeneous distribution of the variables is unstable according to oscillatory 
perturbations inside the wide interval of n 1; which is much wider then interval 
( — 1, 1). At the same time, smooth distributions at the maximum and mini- 
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mum values of n\ are ±^/3 correspondingly. In the first approximation, these 
smooth regions of the dissipative structures resemble homogeneous ones and 
are located inside the instability regions. As a result, the unstable fluctuations 
lead to homogeneous oscillations, and the dissipative structures destroy them- 
selves. We can conclude that oscillatory modes in such type FODEs have a 
much wider attraction region than the corresponding region of the dissipative 
structures. 

It should be noted that the pulsation phenomena of the dissipative structures 
is closely related to the oscillation solutions of the ODE (Fig. 3). Moreover, the 
fractional derivative of the first variable has the most impact on the oscillations 
emergence. It can be obtained by performing a simulation where the first 
variable is a fractional derivative and the second one is an integer. It should 
be emphasized that the distribution of n 2 , within the solution, only shows 
a small deviation from the stationary value (that is why this variable is not 
represented in the figures). 



5 Conclusion 

In this article we consider possible solutions of reaction diffusion system with 
fractional derivatives. Special attention is paid to FODEs linear theory of 
instability which is analyzed in detail. It was shown that three parameters: 
fractional derivative index a the ratio of the characteristic times Ti/r 2 , and 
homogeneous solution n x determine 3-d marginal surface. Inside this surface 
the system is unstable and outside it is stable. Nonlinear dynamics of FODEs 
is investigated by computer simulation of the characteristic examples. 

By the computer simulation of the fractional reaction-diffusion systems we 
provided evidence that pattern formation in the fractional case, at a less 
than a certain value, is practically the same as in the regular case scenario 
a — 1. At a > ao, the kinetics of formation becomes oscillatory. At a = ao, 
the oscillatory mode arises and can lead to nonhomogeneous or homogeneous 
oscillations. 
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